Retrieving Atmospheric Dust Opacity on Mars by 
Imaging Spectroscopy 



S. Doutei 

^Institut de Planetologie et d'Astrophysique de Grenoble (IPAG), France 
(sylvain.doute@obs.ujf-grenoble.fr Phone: +33 4 76 51 41 71 Fax +33 4 76 51 41 46) 

X. Ceamanos^ 
'^Meteo Prance CNRM/GMME/VEGEO 

T. Appere^ 
^Laboratoire AIM CEA-Saclay 



Abstract 

We propose a new method to retrieve the optical depth of Martian aerosols 
(AOD) from OMEGA and CRISM hypcrspectral imagery at a reference wave- 
length of 1 micron. Our method works even if the underlying surface is com- 
pletely made of minerals, corresponding to a low contrast between surface and 
atmospheric dust, while being observed at a fixed geometry. Minimizing the 
effect of the surface reflectance properties on the AOD retrieval is the second 
principal asset of our method. The method is based on the parametrization of 
the radiative coupling between particles and gas determining, with local altime- 
try, acquisition geometry, and the meteorological situation, the absorption band 
depth of gaseous CO2. Because the last three factors can be predicted to some 
extent, we can deflne a new parameter /3 that expresses speciflcally the strength 
of the gas-aerosols coupling while directly depending on the AOD. Combining 
estimations of (3 and top of the atmosphere radiance values extracted from the 
observed spectra within the CO2 gas band at 2 /xm, we evaluate the AOD and 
the surface reflectance by radiative transfer inversion. One should note that 
practically jS can only be estimated for surfaces spectrally dominated by miner- 
als and/or H2O ice. Validation of the proposed method shows that it is reliable 
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if two conditions are fulfilled: (i) the observation conditions provide an airniass 
that is large enough j/ > 3.5 (ii) the aerosol are vertically well mixed in the 
atmosphere. Experiments conducted on OMEGA nadir looking observations as 
well as CRISM EPF produce very satisfactory results. Finally in a companion 
paper the method is applied to monitoring atmospheric dust spring activity at 
high southern latitudes on Mars using OMEGA. 

Keywords: Mars; Atmosphere; Radiative Transfer; Aerosols; OMEGA; 
CRISM 



Introduction 

Visible and near infrared imaging spectroscopy is a key remote sensing tech- 
nique to study and monitor the planet Mars. Although its atmosphere is much 
fainter than Earth's, its composition dominated by CO2 gas implies numerous 
and strong absorption bands that often overlap with spectral features coming 
from the surface. Furthermore small mineral particles or H2O ice clouds of- 
ten drift over Martian surfaces at various altitudes. These aerosols have also 
a strong, spatially and temporally varying influence on the morphology of the 
acquired spectra. As a consequence accurate analysis for the study of surface 
materials requires the modeling and the correction of the atmospheric spectral 
effects. The first step in this matter consists in retrieving the aerosol optical 
depth (AOD) over the scene. 

Since 2004 the imaging spectrometer OMEGA aboard Mars Express per- 
forms nadir-looking and EPF (Emission Phase Function) observations in the 
VIS (visible) and the SWIR (short wave infrared) (920-5100 nm at 14 to 23 
nm spectral resolution) for the study of the surface and the atmosphere of the 
red planet. The spatial resolution of OMEGA typically varies between 300 and 



2000 m/pixel due to its eccentric orbit. The authors in Vincendon et al. (20071 
have developed a method to quantify the contribution of atmospheric dust in 
SWIR spectra obtained by OMEGA regardless of the Martian surface compo- 
sition. Using multi-temporal observations at nadir with significant differences 
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in solar incidence angles, they can infer the AOD and retrieve the surface re- 
flectance spectra free of aerosol contribution. However, this method relies on the 
very restrictive assumption that the atmosphere opacity remains approximately 
constant during the time spanned by the employed acquisitions. 



In Vincendon et al. (20081 the same authors map the AOD in the SWIR 
above the south seasonal cap of Mars from mid-spring to early summer. This 
mapping is based on the assumption that the reflectance in the 2.64 /im satu- 
rated absorption band of the CO2 ice at the surface is mainly due to the light 
scattered by aerosols above most places of the seasonal cap. In this case, one 
geometry is sufficient for the AOD retrieval. Nonetheless, this method is re- 
stricted to the area of CO2 deposits that are not significantly contaminated by 
dust nor water ice. 

The Compact Reconnaissance Imaging Spectrometer for Mars (CRISM) is 
a hyperspectral imager on the Mars Reconnaissance Orbiter (MRO) spacecraft. 
In targeted mode, a gimbaled Optical Sensor Unit (OSU) removes most along- 
track motion and scans a region of interest that is mapped at full spatial and 
spectral resolution (18 or 36 m/pixel, 362-3920 nm at 6.55 nm/channel). In 
the targeted mode, ten additional, spatially binned images (180 m/pixel) are 
taken over the same region before and after the main image at 10 emergence 
angles ranging from -70° to 70°. They provide the so-called Emission Phase 
Function (EPF) for the site of interest that is intended for atmospheric study 
and correction of atmospheric effects. 

Regarding the atmospheric correction of CRISM data, [McGuire and 14| 



co-authors ( 2009 1 adapted and improved the so-called volcano-scan technique 



(Langevin et al. 20061. This method removes the CO2 gas absorption bands 
of any spectrum of interest after division by a scaled reference spectrum (i.e. 
the ratio between the atmospheric transmission at the summit and the base 
of Olympus Mons evaluated on a Martian sol when the amounts of ice and 
dust aerosols were minimal). This simple technique does not correct for aerosol 
effects and cannot operate for icy surfaces. 



In McGuire and 23 co-authors (20081, a DISORT-based model retrieved 
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the dust and ice AOD, the surface pressure and temperature from previous 
experiment products as well as the acquisition geometry, and the measured I/F 
spectrum as inputs. Then, a surface Lambertian albedo spectrum is computed 
as the output. However this algorithm does not retrieve the AOD directly from 
the images nor does it take advantage of the EPF measurements of the CRISM 
targeted mode. 



Brown and Wolff ( 2009 1 proposed a first attempt in that direction by using 



the DISORT algorithm to model the signal at one wavelength (i.e. 0.696 fim). 
They iteratively adjust three parameters (surface albedo, dust and ice opacity) 
in order to achieve a close fit at five points spread across the EPF curve. Never- 
theless the method is time consuming and the surface albedo is assumed to be 
Lambertian. It has been proved that this assumption bias the AOD and surface 



estimation (Lyapustin 19991. 

In this article, we propose an original method that overcomes the previous 
limitations to retrieve the optical depth of the Martian dust from OMEGA or 
CRISM data at the reference wavelength of one micron. The method is based 
on a parametrization of the radiative coupling between aerosol particles and gas 
that determines, based on the local altimetry and the meteorological situation, 
the absorption band depth of gaseous CO2. We consider the intensity of the 
absorption feature at 2 /im as a proxy of the AOD, provided that the other 
influencing factors have been taken into account. When processing OMEGA 



observations, we are complementary to the method of Vincendon et al. (20081 
since our approach processes pixels occupied by a wider variety of materials - 
pure mineral or water ice as well as CO2 and H2O deposits contaminated by a 
large amount of dust - while being observed at a fixed geometry. When process- 
ing CRISM observations we take full advantage of the top of the atmosphere 
spectral EPF measured by the instrument while minimizing the infiuence of the 
surface non Lambertian behavior on the retrieval of the AOD. 

This article is organized as follows. In section [T] we describe the post- 
processing and the formatting of the sequence of EPF calibrated image cubes 
accompanying the high resolution CRISM observation. In section [2] we give 
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insights about Martian atmospheric properties and radiative transfer (RT) in 
the SWIR range. Furthermore we describe models that calculate the spectral 
radiance coming from Mars and reaching the sensor at the top of the atmo- 
sphere. In section [3] we expose the basic assumptions, RT parametrization, and 
properties on which relies the method for retrieving the AOD from the images. 
In addition we describe how we implement the method that is then validated 
on synthetic data. Finally, in section [4] we present and discuss results obtained 
for representative observations, one from OMEGA and the others from CRISM. 
Conclusions are eventually drawn in the last Section. 

1. Post-processing and formatting of CRISM EPF observations 

Due to the composite nature of CRISM EPF acquisitions, such data need 
some post-processing prior to AOD retrieval. First of all, CRISM data are 
transformed from apparent I/F units (i.e. the ratio of reflected radiance to 
incident intensity of sunlight) to reflectance units. A Lambertian surface is 



supposed and data are divided by the cosine of the solar incidence angle (Murchie 



and 49 co-authors[ |2007[ ). 

CRISM data are corrected for artifacts caused by non-uniformities of the 
instrument, residuals of the radiometric correction or external sources. First, 
hyperspectral images are corrected for striping and spiking effects. Correc- 



tions for both distortions have been proposed by Parente ( 2008 1 . Secondly, the 
CRISM spectrometer is affected by a common artifact to "push-broom" sen- 
sors, the so-called "spectral smile" effect. Smile effects refers to the artifacts 
originated by the non-uniformity of the instrument spectral response along the 
across-track dimension, i.e., the horizontal axis that corresponds to the data 
columns. Ceamanos et al. developed a twofold method that corrects CRISM 
observations for smile by mimicking an optimal smile-free spectral response 



(Ceamanos and Doute 20101. In this method, data are first resampled to a 
set of reference wavelengths that correspond to the detectors owning the best 
spectral response. Then, an increase in spectral resolution is performed on data 
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that are particularly affected by spectral smile effects by a spectral sharpen- 
ing strategy. The correction for smile effects is crucial since the spectral band 
corresponding to the absorption maximum of the CO2 gas at 2/im is particu- 
larly affected by spectral smile. We remember that the proposed AOD retrieval 
method is based on this spectral feature. As a matter of fact, the estimation of 
the AOD might be biased if smile effects are not corrected. 

CRISM observations are spatially rearranged to evaluate the atmosphere 
optical depth of a given terrain position depending on geometry. First, the 
central image that owns the highest spatial resolution is binned by a factor of 
ten to match the spatial resolution of the EPF series. After that, all images are 
projected onto the same geographical space. At this point, CRISM observations 
showing a high overlap of the eleven acquisitions can be discriminated. In fact, 
a good overlap assures the existence of a high number of terrain units that 
have been sensed from many geometries. In a single hyperspectral data set 
that gathers all EPF images and the central scan (hereafter called CSP, Cube 
SpectroPhotometric), each super pixel conjugated to a given terrain unit gathers 
up to eleven spectra depending on the completeness of the overlap. More details 



can be found in (Ceamanos et al. 2013). 



2. Mars atmospheric radiative transfer 

In the present section we give insights about Martian atmospheric properties 
and radiative transfer (RT) in the SWIR range. Furthermore we describe models 
that calculate the spectral radiance coming from Mars and reaching the sensor 
at the top of the atmosphere. 

We consider the Martian atmosphere as a vertically stratified medium with 
a plan parallel or spherical geometry. We divide the atmosphere into « 30 
layers spanning the 0-100 km height range. The lower boundary of the layer 
I is at height above the surface. Each layer contains a mixture of gas and 
solid particles and is homogeneous in temperature T; [K], pressure Pi [mbar] 
and gaseous composition C; [ppm]. 
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2.1. Gas 

A dozen atmospheric absorption bands due to gaseous CO2 and marginally 
to H2O are present in the SWIR. We calculate the transmission function of the 
atmosphere along the vertical, for each pixel of our spectral images, and as a 
function of wavelength. 

LBLRTM. For that purpose, we employ a line- by- line radiative transfer model 



(LBLRTM) (Clough et al. 1995) fed by the vertical compositional and ther- 



mal profiles predicted by the European Mars Climate Database (EMCD) For- 



get et al. (20061 for the dates, locations, and altitudes of the observations. 
The LBLRTM code takes into account the vertical profile of temperature, pres- 
sure and gaseous composition and interrogates the spectroscopic database HI- 



TRAN ( «high-resolution transmission molecular absorption database)), http: 



//www ■ cf a . hcirvard . edu/hitrcLn/| ) to compute for each homogeneous layer / 



the monochromatic absorption coefficient aj' [m^] (Ai/ = IG^'^cm^^) of the 
gaseous mixture for the wavenumber v. The absorption lines typical of each 
molecular specie are characterized by a «Voight» profile that describes at the 
same time the broadening by Doppler effect and by collision. The code LBLRTM 
from the Atmospheric and Environmental Research (AER) company also takes 
into account various spectroscopic coupling between species (some of them 
not being directly active in the optical sense like N2 or the inert gases). If 
ui [m"^] is the vertically integrated molecular density in terms of number in 
layer I, the vertical monochromatic transmission follows the Beer-Lambert law: 
T'^(ui) =exp(-a>0 

Method k-correlated coefficients. In order to take into account the fine details of 
the gas lines, the final RT calculation, which will integrate the gaseous as well 
as particulate absorption and scattering, should be conducted in principle ac- 
cording to the maximum spectral resolution of the problem (Ai^ = lO^^cm"^). 
Since we must simulate tens, even hundreds, thousand spectra with numeri- 
cally intensive codes, such a requirement cannot be satisfied practically. Each 
spectral channel k of the sensors that we consider (see Introduction) presents a 
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width « 10 cm"^ and thus contains hundreds of hues that cannot be considered 
individually. Consequently, only a statistical approach like the one proposed by 



Goody at al. (19891 offers a practical solution of the problem. Starting from the 
line spectra in the spectral interval [i^k, i^k+i] {i^k, being the lower bound of the 
current channel k), we build the discrete form of the frequency distribution for 
the absorption: 

rk( ^_ W{a^,ai + Aai) 

J^k+1 - Vk 

Qi , i = 1, . . . , A'^ is a partition of the interval of values taken by the monochro- 
matic absorption coefficient in [vk, i^k+i]i W{ai^ ai + Aa^) the number of spec- 
tral points having a value falling in the bin i of lower bound ai and of width 

First, we choose an irregular partition that gives the same number of points 
for each bin i. The objective is to estimate numerically the frequency f^{ai) (i.e. 
a probability) with bins all containing the same number of individuals (w 50), a 
number large enough to guaranty a correct and homogeneous accuracy. Because 
there are much more lines of weak intensity than lines of strong intensity, the 
width Afli is a strongly increasing function of i. From f^{ai) the discrete form 
of the cumulative frequency distribution g^{a) is derived such as: g^{an) = 
'Y^^^i f^{ai)Aai. The inverse function a(g) = (g*^) ^ (a), which is generally 
defined in the literature as a «k distribution)) and is established in the interval 
g G [0,1], is the highest absorption value reached by the lower fraction g of 
the population of lines ordered according to their intensity. In practice, the 
previous function is calculated for any value of g by spline interpolation of the 
experimental curve passing through the points (5'°(a„),a„). 

Second, we split the line population in {N + 1) sub-sets each representing 
a fraction Ag^ = gi+i — gi , i — I, . . . , N ~ 1, Ago = gi and Agjy = 1 — g^ of 
the total. A Gaussian quadrature gi , i = 1, . . . , N of order N=2™, m e N for 
example allows an optimal sampling of the weak and strong absorption values. 
Then, the function a(g) allows to build the corresponding partition but in the 
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space of absorption values = a{gi). We calculate the mean value of the 
monochromatic absorption coefficient in each bin [ai, ai + Aai] thanks to the 
values returned by LBLRTM AER. 



Finally, as in Lacis and Oinas ( 1991 1, we assume that the mean transmission 



value for the instrument channel k through the layer / is: 



Pl N+l 

T^{ui)= / exp{-aui) f''{a)da ^ / exp {-a{g)ui) dg fa exp(-a-'^u;)Ag^ 

(2) 







Furthermore, we introduce a mean optical depth linked to the sub-sets i = 
0, . . . , + 1, each representing a fraction Ag; of the population of lines sorted 
by increasing intensity : t^-''' = a^^ui. Then T^^^^ = Yli=\ ^» exp {^Tg^^f^ 
with Wi the set of coefficients associated to the Gaussian partition gi. 

Spectral transmittivity. If we now consider the whole stack of atmospheric layers, 
we must make an additional hypothesis in order to calculate an approximation 
of the total vertical transmission. It is based on the fact that the atmospheric 
composition does not change significantly over the heights z spanned by the 
region that contributes predominantly to the transmission. Thus, the form of 
the distribution a{g) is very similar from one layer to the other even though the 
absolute level -which depends on the pressure level for example- varies. Indeed, 
if an intense absorption line is centered at wavenumber v for layer I then the 
line originating from the same transition has a high probability of falling nearly 
at the same position for the layer I' ^ I even though its width and intensity 
are different. In other words, the probability that a strong absorption level 
in the layer I at wave number v is associated to a weak level in the layer I' 
is very low. Statistically such a situation translates to the following rule con- 
cerning the probability that the absorption level for the wavenumber i> belongs 
to the sub-set Ag.; in the layer I and to the sub-group Ag^ in the layer I' : 
P (j^AgiY I {AgjY ^ w Sij (Kronecker symbol). The transmission through both 
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layers that should be generally written: 

N N 
i 3 

(3) 

can then be approximated by : 

N 

T^{uuui,) « exp (-{at-^ui + at'^^uv)^ Ag, (4) 

i 

since we use the same quadratures for both layers. This formula can be gen- 
eralized to the whole stack of layers and the vertical transmission through the 
atmospheric gases in a multilayered medium is given by: 

AT+l 

T^^,,= ^A5,exp(-^4'=_) (5) 

i I 

Because the complete calculation is very time consuming, it is only performed 
for a limited (< 20) number of reference points representative of "regions" sliced 
according to regular bins in latitude lat, longitude long and altitude h. Three 
transmission spectra are generated for the maximum, mean and minimum alti- 
tude of a given region. Then the transmission spectrum of all the pixels belong- 
ing to the region is interpolated from the triplet depending on their individual 
altitude to give Tg^^{h, lat, long) 

2.2. Aerosols 

Optical properties of Martian mineral aerosols are still largely undocumented 



even though recent studies have improved our understanding (Korablev et al 



20051. We favor the single scattering albedo, optical depth spectral shape and 



phase function retrieved in the near-infrared by Vincendon et al. ( 2008 1 and 



Wolff et al. (20091. In the first case, the properties correspond to a specific 



population of aerosols taking place at high latitudes in spring, our preferential 



field of investigation (see section 4.2 1. Nevertheless the phase function -a coarse 
Henyey-Greensten function with one parameter- is valid only for the phase angle 
range spanned by the data set of nadir OMEGA observations that we consider. 
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In the second case, since it has been proved that the Henyey-Greenstein function 
is somewhat inaccurate to recreate the effect of aerosols for the CRISM angu- 
lar ranges, the radiative properties are described with more details. A single 
scattering albedo, optical depth spectral shape, and phase function based on 
cylindrical particles are considered. Nonetheless such radiative properties are 
not specific of the polar regions. The nature of the mineral particles drifting in 
the atmosphere should not change fundamentally with latitude, but their de- 
gree of hydration certainly does . Fortunately, for the typical hydration values 
reached on Mars, the spectral effects are only noticeable between 2.7 and 3.5 fim 



(Pommerol et al. 20091, a range that we thus ignore. Water ice crystals may 
also be frequently present but are not considered in this paper because they may 
influence the CO2 gas 2/im absorption feature in a specific manner. Because the 
atmosphere is usually well mixed, we take an exponentially decreasing number 
of aerosol particles of constant characteristics with height such that the vertical 
distribution of optical depth at a reference wavelength (fco : 1 /^m) is given by : 

ra'er° = T^^A^M^z' / H scale) - exp{-z' + ' / H scale)) 
/(I - CXp{-Zmax/H scale)) 

'Taero = X^^^i '^aer '^^ ^hc column integrated opacity of the aerosols at the refer- 
ence channel /cq, H scale is the scale height of the distribution, z' the height of 
layer I lower interface and zmax — 100 km. The scale height Hscaie and the 
optical depth t^^^ remain the only free parameters concerning the atmosphere. 

2.3. Gas-aerosols coupling 

In general, photons undergo absorption and multiple scattering events when 
interacting both with gas and aerosols. One should note that we neglect Rayleigh 
scattering due to atmospheric molecules because, according to the formula of 



Bodhaine et al. (19991, that will lead for a mean ground pressure of 6 mbar 



to an optical depth of Tray ~ 5.10^'^ - far lower than the contribution we can 
expect for the aerosols in any case. The value of the monochromatic aerosol 
optical depth (at Av = 10^*cm~^) does not vary significantly within a channel 
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Gas 


Aerosols 


Single scattering albedo 


^gas — 


^aer 


Optical depth 


Tgasi ,^gi, i = 0, . . . ,7V + 1 


^aer 


Phase function 


T — n 

^ gas ^ 


T 

aer 


Tt - Tgas i + Taer ; '^i - ^aer 1 (1 + ) , I - 0, . . . , iV + 1 T - T aer 



Table 1; Sum up of the physical quantities characterizing the "elementary" processes of the ra- 
diative transfer taking place in each homogeneous atmospheric layer. We omit in the notation 
the indices I and k for readability. 



k whereas the counterpart for the gas can vary by several orders of magnitude. 
As a consequence calculating at channel k the transfer of radiation in the atmo- 
sphere requires to couple gas and aerosols properties at each atmospheric level 
/ and for each sub-set i = 0, . . . , iV -f 1 of the gas line population. As indicated 
in Table [l] the coupling procedure leads to -I- 1 sets of optical depth, single 
scattering albedo and phase function that characterize the layers. These prop- 
erties, as well as the acquisition geometry (9i, 9e, (f'e), are introduced one by one 
in a RT engine (e.g. DISORT, Stamnes et al. ( 1988[ )) in order to calculate for 



each level z the partial radiance field Oi, 9e, 4>e) corresponding to each sub- 
set lines+aerosols. The surface is characterized by its bidirectional reflectance 
distribution function (BRDF). The total radiance field is the mean value of the 
partial fields weighted by the fractions Agi : 

N+l 



3. Method for retrieving the optical depth 

The proposed method is based on a parametrization of the radiative coupling 
between aerosol particles and gas that determines, with local altimetry and 
the meteorological situation, the absorption band depth of gaseous CO2. The 
coupling depends on (i) the acquisition geometry (ii) the type, abundance and 
vertical distribution of particles (iii) the bidirectional refiectance factor of the 
surface (BRF). For each spectrum of an image, we compare the depth of the 
2 /im absorption band of gaseous CO2 that we estimate on the one hand from 
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the observed spectrum and on the other hand from a calculated transmission 
spectrum through the atmospheric gases alone. This leads to a relevant new 
parameter that directly depends on r^g^. Combining the latter with the radiance 
factor R^lg within the 2 /^m band, we evaluate by RT inversion the AOD t^^^ 
and the reflectance factor of the surface. 

3.1. Parametrization of the spectral signal 

Our method is based on two main assumptions: 

1. First we assume that the surface reflects the solar and atmospheric ra- 
diation isotropically and. consequently, that the ground-atmosphere interface is 
characterized by a single scalar quantity (the normal albedo A'^^^j) that depends 
on wavelength. 

2. Second we parametrize the radiative coupling between aerosols and gas 
by assuming that the latter contribute to the signal as a simple multiplicative 
filter. In this way, the top of the atmosphere (TOA) radiance is written such 
that 



l''{e,,9,,^,)^T^^,{hJatJongY''<^'^-'-'f'-^^^^^^^^^ 

(7) 

the effect of the acquisition geometry, the aerosols, and the surface Lamber- 
tian albedo A^^^y being taken into account by scaling the aerosol free vertical 

gaz V 

at: 

surf+aer V 



transmission T^^^(h,lat,long) -calculated according to 2.1- by an exponent e. 



^surf+aeri^i^ 'Pe) IS the Calculated TOA radiance but without considering the 



atmospheric gases in the radiative transfer. Factor e'" can be further decomposed 
into two terms: 



with the airmass i/ ~ cos{0 ) ^ cos(6> ) ■ F^^ctor tp is purely geometric and allows 
a quick and simplified calculation of the free gaseous transmission along the 
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Po 


Pi 


P2 


P3 


P4 


OMEGA 


1.2 


0.66312127 


0.44429186 


-0.024559039 


0.00068174262 


CRISM 


0.95 


0.66312127 


0.44429186 


-0.024559039 


0.00068174262 



Table 2: Value of the coefficients allowing to calculate 



acquisition pathlength (incident-emergent) using the vertical free transmission 



= f: a. exp(- fit ..^) « {Tl.)^^"^ (9) 

i I 

the approximate formula being = Po{pi+P2t^+P3t^'^ +Pii^^)- The values 
of the previous coefficients have been calculated by fitting the second term of 
equation |9]by the third term and are given in Table [2] Factor po accounts for the 
difference of spectral resolution and radiometric calibration between OMEGA 
and CRISM. 

Factor /3'^' on the other hand expresses the aerosol effect on gaseous ab- 
sorption and depends on the considered channel k. In practice this spectral 
dependability can be easily managed. Indeed, numerical experiments and the 
analysis of real data (see following sections) show that this exponent depends 
principally on the gaseous absorption intensity Tg^^. If (i) /3'' can be estimated 
experimentally (13'^) for one or for several geometries for each pixel (super-pixel) 
of an OMEGA (CRISM) observation and (ii) if one assumes a value for A^'^^.^ 
(fci is the spectral band number at 2/im) then {rj^^r) '^^'^ be derived. Indeed, 
function l3''{0i,0e,(t>e,Tl^erT ^scaie, A'^u^j) is invcrtiblc provided that the scale 
height Hscaie of the aerosols is known. In this matter, several studies such as 



Vincendon et al. (20081 suggest that Hscaie = 11 km is a good guess, but in 



section [4] we prove that Hscaie = 8 km may be more appropriate as it leads to 
better results. 

3.2. Numerical aspects LUT 

Our optical depth retrieval system uses a multi-dimensional look-up table 
(LUT) of j3 values arranged according to discrete combinations of acquisition 
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geometries {6i, Oe,(t)e) and of physical parameters (r^g^, H scale, ^s^/)' Table 
[Sj For each combination, the value of /3 is calculated on the basis of three TOA 
spectra. T^^f ^ reference free vertical transmission corresponding to a given 
time and location on the surface of Mars. It has been chosen in order to cover a 
range of absorption values (all spectral bands considered) typical of the southern 
hemisphere in spring. The spectra J*^ and Igurf+aer units of radiance are 
generated successively by the code DISORT with and without considering the 
gaseous absorption in the atmosphere. A series of k-correlated coefficients f^^*^ ^ 
intervenes in the calculation of T^'^j, ^ and l'^ . Finally the LUT is built according 
to the formula: 

k=k[ k=k[ surf+aer 

(10) 

with fk — l/K. The channels k = k'l, . . . ,k'{ {K = k'{ — k'l + 1) encompass 
the 2 /im CO2 gas absorption band such as T^^^ < 0.85. For DISORT calcula- 
tions, we consider the atmosphere as a plan parallel vertically stratified medium 
which is a sufficient approximation to calculate the TOA radiance whenever 
incidence and emergence directions are «10° above the horizon. Several exper- 
imental scatter-plots T^^^ - (5'' built using the reference and also other sets of 
k-correlated coefficients show that there exists a linear relationship between the 
two quantities with an excellent correlation factor (> 0.9 see figure [T| : 

Z?*"' = a{l', T^g^, H scale, A^su^f) + ^(^J '''aer, Hscale, A'lmr f)T^ref (H) 

We notice that the coefficients (a, h) of the regression depend weakly on the 
parameters T^°^,Hscaie,A^J^^f. 

3.3. Estimation of the coupling factor gas-aerosols 

Because the intensity of the CO2 gas 2/im absorption feature is diversely 
accessible from the spectra depending on the nature of the surface materials 

^the albedo of the surfaces that we consider rarely exceeds 0.6 
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Table 3: Series of values taken by the incident, emergent and azimuth angles as well as by the 
physical parameters used for the construction of the reference LUT 
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Figure 1: Experimental scatter-plots between T^^j^ and f}^ for different geometries and atmo- 
spheric conditions. The magenta and cyan points, fitted by regression lines of the same color, 
come from synthetic data. The blue, green and red points fitted by a single regression line (in 
red) come from the estimation of factor Z?'' for different observed spectra (see section 3.3 I 



16 



present in the pixel, we need to consider different strategies for estimating the 
factor f3. 



Mineral surfaces: method 1. The factor /? can be readily estimated for every 
OMEGA or CRISM spectrum that does not show any H2O or CO2 ice signature 
superimposed on the 2 /im CO2 gas absorption band. In this way, this absorption 
band is completely due to the atmosphere. A similar formula to the one used in 



the volcano scan technique ( McGuire and 14 co-authors 2009 1 is used replacing 



the Olympus reference transmission spectrum by (Tg^^^)'^^'^^ such that 

pt" = 



ipiy) In 



.gas 

'-gaz 



In ( |fM + 0.0909 In 



(12) 



with fca the index of a channel falling at the shortwave extremity of the 2 
/im band wing. The scatter-plot T^^^-fS'', k = k[, . . . ,k'( built for each pixel 
of a real OMEGA observation confirms the existence of a linear relationship 
between the two quantities with a generally moderate dispersion (see Figure 
[1]). The regression coefficients (a,S) are then calculated and will allow us to 
estimate exponent /3 in the same absorption range than the reference ; 

/3,e/ = ^ + ^ E A-^-/ (13) 

Note that we use the series of channels k = k[, . . . ,k'( for the estimation 
of Pref allowing us to be less sensitive to noise and channel aging than using 
the single channel ki. In order to perform the inversion process, we take as an 

k R^^ 

initial value for the lambertian albedo of the surface A ^ — "'" 



surf (Tgi^)^M^''^ ' 

Water ice surfaces: method 2. A specific procedure is necessary for any spec- 
trum marked by the signature of water ice with the possible presence of dust 
but without spectral contamination by CO2 ice. A numerical optimization is 
performed regarding a cost function that depends on /3 and that expresses the 
quality of the gaseous absorption correction on the spectra. The quality criteria 
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we choose is the band shape of solid H2O at 2 /im that must show a unique local 
minimum and a simple curvature. Mathematically this criteria is equivalent to 
the second derivative of the spectrum that must be close to zero for channels 
[k'l, k'{] covering the H2O band and slightly beyond. After all the cost function 
is written: 

^^^^^£;(^(^|rf^) ^^^^ 
Searching for the global minimum of this function leads to the estimation of 
the exponent /3 : l3 — argmin^(/3). The normalization of /3 up to the reference 
level of absorption uses the regression line: 

Once again, in order to realize the inversion we take as initial value for the 
lambertian albedo A^^^^ = (t^Hm^ ' 

Carbon dioxide ice ( C02)- Currently there is no simple strategy to estimate in a 
pixelwise manner the /3 coefficient from spectra presenting the signature of solid 
CO2 because the latter cannot be distinguished from its gaseous counterpart at 
2/im. For this kind of surface, our method to estimate r^^^ is ineffective. 

If, on the one hand, CO2 ice proves to be sufficiently pure with grain sizes 
more than 100 /im in diameter (this conditions is always respected in practice 



for the seasonal deposits) then the method proposed by ( Vincendon et al. 2008 1 
is the only solution. It is based on the assumption that the reflectance level at 
2.63 /im (channel fc2) is nearly equal to zero because of the absorption saturation 
of the ice. Consequently, in this case, the reflectance factor R^l^ represents the 
path radiance of the atmosphere (mostly due to the aerosols at channel ^2) 



and is a bijective, easily invertible function of t^^^. In section 3.5 we will see 
how exponent /3 can be calculated from t^^^, computed using (Vincendon et al. 



2008 1 , and R^l^ for correcting the gaseous absorption on the spectra (section 
3^51. 
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If, on the other hand, CO2 ice proves to be substantially contaminated by 
dust or water ice, then the only conceivable way to estimate (3 could be a 



statistical procedure as it is proposed in Luo et al. (20101. 

Hereafter, we simplify the writing of the parameters by removing subscript 
ref from the notations. The quantities /3 and /3 are as if they were wavelength 
independent because they are normalized up to the reference level of absorption. 

3-4-. Sensitivity and feasibility study. 

Synthetic data. We now examine the sensitivity of function /3 regarding its 
different parameters di,6e, (f>e,'''aerT -^scaie^ ^surf- This study is summarized as 
follows. The exponent /3 decreases with the airmass v especially since aerosol 
opacity is high. This is the direct consequence of an intensifying coupling ef- 
fect between gas and aerosols, the latter reducing the effective pathlength of 
the photons in the atmosphere by scattering and thus the probability to be 
absorbed by CO2 gas. The exponent evolution reflects such a decrease of rela- 
tive absorption band intensity. The coupling is higher when the sensor looks in 
the solar direction (90 (f>e 180°) than when its looks in the opposite half 
space (0 < (/>e ^ 90°) . The higher the scale height characterizing the vertical 
distribution of the aerosols, the more aerosols influence atmospheric absorption 
especially since their opacity is high. Finally the surface reflectance factor also 
controls the absorption by CO2 gas because of occurrence of multiple scatter- 
ing between the two media. A bright surface promotes a higher number of 
roundtrips for the photons which then have a higher probability to be absorbed: 
exponent f3 increases with Lambertian albedo A'^l^^j:. 

The experimental behavior of l3, evaluated above mineral surfaces based 
on real data with respect to airmass i' and r^g^, is also full of information. 
We distinguish two situations according to whether we work with OMEGA or 
CRISM observations. 

OMEGA data. In the OMEGA case, we consider large populations of pixels 
extracted from several hyperspectral images. The principal variability of the 
aerosol opacity across the scene comes from altitude changes -sometimes several 
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kilometers in elevation- that modify the local atmospheric height and thus the 
column integrated density of the aerosols. The horizontal heterogeneity of the 
particle density at constant altitude generally implies a more moderate vari- 
ability. Hence, by selecting pixels belonging to slices spanning no more than a 
couple hundred meters in elevation but observed under various geometries, we 
can draw a scatter-plots {i/, (3) such as that presented in Fig. [2] Set apart a no- 
ticeable internal dispersion linked to factors such as those previously mentioned 
or spatial variations of Lambertian surface albedo, the points corresponding 
to the OMEGA observation 1880_1 are organized according to monotonously 
decreasing curves of /? over a wide range of increasing v values. The curves 
faithfully follow the theoretical lines that can be built thanks to the numerical 
reference LUT for H^caie w 8 km and A'^l^j.- w 0.3. The best fit gives use m 
each case an estimation of the mean optical depth r^^^ . However we notice that 
the experimental curves display a systematic disagreement with the described 
theory for 1/ < 3.5 because of a sudden increase of f3. 
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CRISM data. The sensitivity study is conducted on CRISM targeted image 
FRT82EB. This observation corresponds to an ice-free cratered area of the south 
pole of Mars (^ai=-83°). The image geometry corresponds to the conditions of 
the high latitudes with 9e=74.5°. Due to the sun-synchronous orbit of CRISM, 
there are two modes in azimuthal angle, ipe —34° and 0e ~146°. Regarding the 
topography, FRT82EB presents mild slopes and a constant altitude. The initial 
albedo of the surface can be approximated by the average value of the spectral 
band at 1 micron. In this case, the A^^l^j of FRT82EB is equal to 0.3. 

Data are processed by the data pipeline described in Section [l] After pro- 
cessing, spectral bands of interest {e.g. IR bands 137 and 154) present a high 
signal to noise ratio and are artifact-free. Thanks to the good overlap of the 
high-resolution image and of the EPF, we can define up to 850 super pixels that 
are represented by more than 6 geometries. 

The factor /3 can be readily estimated for every CRISM spectrum that does 
not show CO2 ice signature superimposed on the 2 jim CO2 gas absorption 
band. Once a series of spectra corresponding to a super pixel is treated, we 
obtain /3 as a function of up to eleven geometries. Figure [3] summarizes the 
results for the image FRT82EB by plotting /3 as the function of the airmass 
distinctively for the two azimuths explored by the observation (black and red 
crosses) . 

One can see how the azimuthal angle has an impact on the gas-aerosol cou- 



pling due to aerosol phase function strong anisotropy ( Wolff et al. 2009 1 . In 
addition, the model curves of /3 that provide the best matching with the experi- 
mental results are shown by plain lines. In fact, the best match corresponds to a 
value of r^g^ around to 0.5 and Hgcaie of 8 km. On other targeted observations, 
contrary to the OMEGA case, we do not observe this time the disagreement for 



u < 3.5 (Fig. 10 1. We do not have a convincing interpretation of this difference 
for the moment. 

The practical conclusion of these experiments is that our method when ap- 
plied to the OMEGA data set is only effective if the airmass of the acquisition 
is large enough (i^ > 3.5), this constrain being partly lightened for the CRISM 
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Figure 3: Experimental /3 curves from the CRISM image FRT82EB according to airmass. 
The model curves that provide the best match are also plotted as solid lines. 



images. In both cases, because the sensitivity of /3 with r^^^ (^) increases 
with the airmass, we can expect significant estimation errors as soon as < 3. 



Indeed the ground pressure provided by the Martian climate database (Forget 



et al. , 2006 ) is only indicative of the real pressure at the moment of the obser- 



vation. We must stress that the ground pressure determines T^^^ then /3 and, 
for this reason, an uncertainty Apsurf will propagate as an uncertainty A/3 on 

dp' 



the coupling factor, then as an uncertainty on the optical depth Ar^^,. = ^A.(3 



. We quantify this effect in the validation section [X6| 
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3.5. Algorithms for the simultaneous retrieval of surface albedo and aerosol 
opacity: 

We now describe the practical implementation of the method that allows the 
computation of surface albedo and aerosol opacity maps from a given observa- 
tion. The implementation comes in two flavors depending on the origin of the 
data. 

OMEGA . The processing of a nadir OMEGA image is conducted on a pixel wise 
basis after segmenting the image according to an automatic detection of the 



surface components Schmidt at al. (2007). For the pixels that do not present 



the signatures of solid CO2 the following iterative algorithm is carried out: 
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Input: a spectrum i?^^^ of rank p in the flattened image, the corresponding 
vertical transmission Tg^^ calculated as explained in section 2.1 the acquisition 
geometry (0^, 9e,4>e) and a global estimation of the scale height H scale 

1. estimating the exponent /3 



(a) using method 1 (section 3.3 1 for a mineral surface 

(b) using method 2 for a water ice surface without CO2 ice contamination 

DO 

building a ID curve r^^^ as a function of ji by quadrilinear interpolation 
of the LUT described in section 



3.3 



for a given set {0i,9, 



A 



fci 



surf 



and H. 



scale ■ 



linear interpolation of the previous curve at /3 to calculate {Taer)i+i 
building a ID curve A^^^j: as a function of R^qj^ by quadrilinear inter- 



polation of the LUT for a given set of 

linear interpolation of the previous curve at 
-late {At^^)^^^ 



) and r^g^ — (r^e^) 

r':l 



+ 1 

in order to cal- 



8. i = i + l 

9. UNTIL 
imax 



^surf 



^surf 



< 0.01 AND 



V aer. 



< 0.01 ORi > 



The algorithm usually converges after a dozen iterations and leads to a complete 
solution once , is calculated from R^i^ and r^iL based on a dedicated LUT. 

S LL f J (JOS LLCt 

If the initial estimation /3 or, to a lesser extent, the estimation of^yl^J^j,^^ is 
poor then the convergence is not reached and there is no solution. 

For the pixels presenting a signature of solid CO2 pure enough (according 



to the criteria established by Vincendon at al. (2008)) we have seen in section 
3.3 that it is possible to obtain rj^^^ by their independent method. The problem 
is then laid down differently: we must calculate the exponent f3 from r^^^ and 



^obs tli^t ^ill used for correcting the gaseous absorption on the spectra. The 



algorithm is also iterative in this case : 
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Input : a spectrum -R^bs '-'^ rank p in the flattened image, the corresponding 
vertical transmission T^az extracted from the ATM cube at the same rank, 



the acquisition geometry 

H scale 



^2 ; ) v^e 



) and a global estimation of the scale height 



1. estimation of r^^^ using the method by Vincendon at al. (20081 

2. initialization i ~ 0, Lambertian albedo ^A^^,,^^ 



0.01 



3. DO 

4. building a ID curve /3 as a function of r^^^ by quadrilinear interpolation 

for a given set 



9. 



of the LUT described in section 
and H scale- 



3.3 



linear interpolation of the previous curve at t^^^ to calculate 



calculation of CL„ = — ; rrfrr^ — 



building a ID curve ^^^^^ as a function R^qa quadrilinear interpo- 
lation of the LUT for a given set {9i, 9e, 4>) and r^^^ 

linear interpolation of the previous curve at Cqj^j, in order to calculate 



surf 



i = i + l 



10. UNTIL 
imax 



surf 



< 0.01 AND 



/3 - P 



< 0.01 ORi > 



CRISM. The processing is applied to the single hyperspectral data set that 
gathers the central scan and the EPF images, i.e. the CSP cube. We choose 
an area of the image that offers the best compromise between its size and its 
spectral homogeneity. In addition, this area should not present the signatures of 
solid CO2. AH the super-pixels that fall in the area are blended into one vector 
of parameters that needs to be explained as a whole by the model provided the 
solution for parameters r^^^ and R^l^ is found. The following iterative algorithm 
is carried out. 
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Input : the collection of valid spectra R^^^ , the corresponding vertical transmis- 
sions Tg^^ extracted from the ATM cube at the same coordinates, the acquisition 
geometries {0i,9e,(t>e) and a global estimation of the scale height H scale 

1. for each valid pixel p, estimating the exponent /? 



(a) using method 1 (section 3.3 1 for a mineral surface 

(b) using method 2 for a water ice surface without CO2 ice contamination 

2. initialization i ~ Q, mean Lambertian albedo for the selected area| A'^^ - ) 

3. DO 

4. building a matrix of /? values, each element corresponding to the 
coupling factor that is predicted by the model for the spectrum of index p 
assuming that the atmospheric opacity is(T^er)^ ,9 = ■ ■ ■ iQ, a sampling 
of the range of possible variation for this parameter. The value is 
calculated by quadrilinear interpolation of the LUT described in section 

for a given set (0., </>e), A^i,/ = (^),' ^^^er) and H,. 



3.3 



scale- 
kO 



5. finding the global minimum of the objective function x{'''aer) — 
Ep=i,...,P - /3^')'to estimate {r^^r),+, 

6. for each valid pixel p, building a ID curve ^^^^^ as a function of i^^A 
by quadrilinear interpolation of the LUT for a given set of (9i, 9e,(t>e) and 
^aer — {'^aer) i+i- Linear interpolation of the previous curve at ^^ti 

in order to calculate ( A^J, 

7. calculation of the mean Lambertian albedo ( A^J^^j ) for the selected 



4 + 1 



area 

8. i^i + l 

9. UNTIL 
imax 



^surf 



< 0.01 ORi > 



3. 6. Validation 

The validation is realized by inverting synthetic data, that is to say, real- 



istic TOA spectra / calculated with DISORT (section 2.3 1 to which we add 
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zero-mean Gaussian noise simulated using a meaningful covariance matrix, i.e. 
calculated from an estimation of OMEGA noise. In the reference simulation, we 
set the properties of the surface A^^^j — 0.3 (minerals) and of the atmosphere 
'''aer — 0-6, H scale = H km, and initial pressure Psurf- In a second simula- 
tion, for the calculation of the spectrum l'', we apply a pressure deviation of 
^Psurf = ±15 Pa typical of the Martian meteorological variability (barocline 
activity, (Forget et al. 1999[ )) but the factor j3 is estimated with a transmission 
spectrum T^az calculated according to the initial pressure. In a third simu- 
lation, we perturb the initial distribution profile of the aerosols (exponential) 
by increasing the opacity of one given layer I, I — 1, . . . in succession by an 
additional amount At^^^ = 0.1. Thus, the total AOD becomes 0.7 with the 
perturbed profile and there are as many runs as the layer number. In each run, 
the geometry varies according to the values given in Table [3] In agreement with 
the sensitivity study (section 3.4 1, Figure [4] shows that our method estimates 
'^aer with an accuracy comparable to the usual methods in the reference case 
provided that the airmass exceeds w 3. Since we do not have access to the 
real ground pressure but to a value predicted by a MGCM, we expose ourselves 
to potentially important errors if < 4 but lower than 10% beyond (4 < < 8). 
Figure [5] illustrates the fact that a perturbation of opacity does not have the 
same impact on factor /3 and thus on the estimation r^^^ according to its height 
above the ground. A low lying aerosol layer is not detectable by our method 
because it does not have much influence on the mean pathlength of the photons 
into the gas. By contrast, a detached layer high in the atmosphere superimposed 
on the standard proflle has a strong influence on the coupling. Thus, such a 
layer could skew our estimation because the latter does not take into account 
such a local perturbation. In conclusion, we state that our method is reliable if 
two conditions are full flUed: (i) the observation conditions provide an airmass 
that is large enough, i.e. v > 3.5 (ii) the aerosol particles are well vertically 
mixed caused by a vigorous convection. 
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Figure 5: Influence on the estimation t^^^ of perturbing the reference vertical distribution 
profile of the aerosols by a local perturbation of opacity 0.1 as a function of the airmass v. 
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4. Results 
4.1. OMEGA 

Figure [6] shows in a synthetic way all the products obtained by the chain 
of operations described above for the observation OMEGA 0RB1880_1 that is 
representative of many. The chosen image covers a large geographical area of 
Mars at high and medium southern latitudes in spring. The bright extended 
part appearing in the maps R^l^ and ^4^°^^ represents the seasonal deposits of 
frozen CO2 as well as the permanent cap that is, at that time, buried by the 
latter. The black coded area indicate the surfaces (among then a significant 
fraction of the "cryptic" region - where CO2 very much contaminated by dust 
exists - for which no method for estimating t^^^ is currently available or indicate 
the pixels for which the iterative inversion has not worked (for example near the 
limb). 

A fast cross validation of the two methods used in order to produce the map 
'''aer ~ o^^'^ ^^^^ and the one by Vincendon et al. (20081- can be performed as 
follow. We plot along-track profiles of optical depth through the maps and we 
examine to which extent the segments corresponding to the mineral surfaces (in 
red) and those corresponding to the frozen surfaces (in blue) are well connected 
in figure |7] The right-side red segment with high dispersion corresponds to the 
"cryptic" region. The quality of the connection seems quite good for 0RB1880_1 
and for other observations that we treated. In addition, by drawing r^^^ as a 
function of the altitude h (figure [8| one can verify the coherence between the 
scale height deduced from the analysis of the observation (red regression by an 
exponential) and the value that we use a priori in the inversion process. 

Finally an enlightening comparison is put forward between the AOD re- 
trieved from our analysis of the near infrared 1880_1 image and an RGB com- 
position made by extracting the TOA martian reflectivity measured by OMEGA 
at three wavelengths (0.7070, 0.5508, and 0.4760 //m) from the corresponding 
visible image (Figure j9]). The composition is stretched so as to make visible 
the dust clouds as yellowish hues against mineral surfaces while the seasonal 
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Figure 7: A MEX along track-profile of r^^^ extracted from the aerosol optical depth map 
1880_1. The method presented in this paper (factor /3) produces the red segment while using 
the method by [Vincendon et ah] l |2008[ ) leads to the blue segment 
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Figure 8: Scatter plot of t^^,, as a function of altitude h for the OMEGA observation 1880_1 
(blue crosses) . The red line is the result of applying a regression to the data with an exponential 
and it corresponds to J/scaie =8 km. 
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Figure 9: Global OMEGA observation 1880_1. The left column displays an RGB composition 
of the TOA martian reflectivity in the visible which is stretched so as to reveal dust in the 
atmosphere as yellowish hues. The central column displays the Aerosol Optical Depth map 
at 1 micron. 



deposits appear completed saturated. Excellent qualitative agreement can be 
noted between the two products down to the details. In the companion paper 
other successful examples are provided. 

4.2. CRISM 

In the CRISM case, six observations of the landing site of the Mars Explo- 
ration Rover "Spirit" (i.e. West of Columbia Hills), acquired with varied ge- 
ometry and atmospheric conditions (see Table |4]) , have been chosen to test our 
method (hereafter identified as the "/3 method"). For that purpose, these scenes 
are of special interest since the AOD values were also retrieved from the same 
data by Wolff et. al. (personal communication) and from PanCam synchronous 



measurements at the ground looking towards zenith ( Johnson et al. 2006 1 . The 
previous procedures are identified as the "W method" and the "P method" re- 
spectively. Special attention is paid by focusing on areas of the observations with 
altitude variations not exceeding one hundred meters and with negligible slope 
at the scale of the pixel size and above. Following the procedure described in 
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Observation FRT 


334D 


7D6C 


812F 


95B8 


B6B5 


3192 


Incidence angle 


55.4 


36.4 


32.6 


39.3 


56.4 


60.4 


Aziniutal modes 


64, 118° 


68, 104° 


88, 97° 


71, 113° 


51, 129° 


65, 135° 
















r^g"^ "W method" (1 ^m) 


0.416 


1.74 


2.07 


0.6 


0.43 


0.38 


RMSE w/model 


1.87E-04 


1.40E-04 


NaN 


1.24E-04 


2.80E-04 


1.39E-03 
















T^°^ "P method" (1 ^m) 


0.871 


1.32 


0.93 


0.66 


0.48 


0.31 


RMSE w/model 


NaN 


1.06E-04 


2.78E-04 


1.44E-04 


2.20E-04 


1.55E-03 
















T^l^^ "/3 method" (1 ^m) 


0.46 


1.41 


1.00 


0.36 


0.46 


0.38 


RMSE w/model 


1.92E-04 


1.14E-04 


2.71E-04 


1.09E-04 


2.41E-04 


1.39E-03 



Table 4: Optical depth at 1 fim retrieved by three difFerent methods for our selection of obser- 
vations. The Root Mean Square Error (RMSE) refers to the adequacy of a non Lambertian 
surface-atmosphere radiative transfer model fed by the latter optical depth in reproducing the 
TOA reflectance factors measured by CRISM. 



section [T] the eleven CRISM images of each FRT observation are combined and 
binned at about 300 meters per pixel. Then, the algorithm adapted to CRISM 
for the simultaneous retrieval of surface albedo and aerosol opacity is applied on 
the CSP cube (section [ij. The operation leads to the results tabulated in Table 
[4] along with those obtained by the "W" and "P" methods. In addition Figure [TO] 
illustrates to which extent our best model matches the experimental (3 versus 
airmass curves for each observation. Similarly to the OMEGA case, airmass 
« 3 represents a threshold above which the fit becomes very satisfactory with 
the exception of the observation FRT95B8. One should also note that the larger 
the difl:erence between the two azimuthal modes is, the more constrained the 
best solution likely becomes because the two branches of the /3 versus airmass 
curves are increasingly separated. 

Once the content of aerosols is known, CRISM TOA radiances can be cor- 
rected in order to retrieve surface BRF. The correction of remotely sensed im- 
ages for atmospheric effects is however not straightforward due to the anisotropic 
scattering properties of the atmospheric aerosols and the materials at the sur- 
face. Traditional inversion algorithms are based on reductionist hypothesis that 



assumes that the surface is lambertian (McGuire and 23 co-authors 2008). Al- 
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Figure 10: Experimental /3 curves from CRISM selected images as a function of airmass plotted 
distinctively for tlie two azimutfis explored by the observation (black and red crosses). In each 
case, the model curves that provide the best match are also plotted as solid lines. 

34 



though this assumption largely simplifies the inverse problem, it critically cor- 
rupts the angular shape of the retrieved BRF since solid surfaces are hardly 
isotropic. Recently, we have proposed an original inversion method to over- 
come these limitations when treating CRISM multi-angle observations |Cea-| 



manos et al. (2013). This inversion algorithm is based on a TOA radiance 



model that depends on the Green's function of the atmosphere and a semi- 
analytical expression of the surface BRF. In this way, robust and fast inversions 
of the model on CRISM TOA radiance curves are performed accounting for the 
anisotropy of the aerosols and the surface. The root mean square error achieved 
by the model when reproducing the CRISM TOA reflectance factor curves is 
indicated in Table |4] for each observation and for each AOD input. A NaN value 
indicates that the inversion was not successful, i.e. no valid surface BRF could 
be retrieved with the proposed aerosol properties and atmospheric opacity. 

By examining Table [4] we conclude that our method gives results generally 
in good agreement with the values retrieved by Wolff et al. with the notable 
exception of observations 812F and 95B8. In the first case, the "W" method 
overestimates r^^^ since, with such an opacity, the inversion algorithm fails to 
produce a physically meaningful surface reflectance. At contrary our value is 
in agreement with the one estimated from the PanCam measurements. In the 
second case, the /3 versus airmass curves are not well fitted by our model imply- 
ing that the estimated value for r^^,. by the "/? method "may not be accurate. 
One should also note that the root mean square error achieved by the TOA 
reflectance factor modeling using our values is lower than when using the values 
given by the "W method" except for observation 334D. Because the coupling fac- 
tor /3 is sensitive to the surface properties through multiple scattering between 
the latter interface and the atmosphere, the retrieval of r^^^ is likely weakly 
dependent on the Lambertian assumption. By contrast, methods that directly 
invert a RT model on the TOA radiance in the continuum of the spectra, such 
as the "W method" are very sensitive to the surface bidirectionnal factor. As a 
result, untangling between surface and aerosol phase dependency is tricky and 



leads to bias if the Lambertian assumption is used (Lyapustin 19991. We ex- 
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plain the moderate discrepancies of results that can be observed between the "P 
method" and the "/? method" if we consider that the first one is more sensitive 
to the low lying layers of aerosols whereas the second one is more sensitive to 



aerosols layers situated at « 2 km in height and above (see section 3.6 1. 
Conclusions 

In this article we propose a method to retrieve the optical depth t^^^ of Mar- 
tian aerosols from OMEGA and CRISM hyperspectral imagery. The method 
is based on parametrization of the radiative coupling between particles and gas 
determining, with local altimetry, acquisition geometry, and the meteorological 
situation, the absorption band depth of gaseous CO2. The coupling depends on 
(i) the acquisition geometry (ii) the type, abundance and vertical distribution of 
particles, and (iii) the surface albedo A^^^.^. For each spectrum of an image, we 
compare the depth of the 2 fim absorption band of gaseous CO2 between (i) the 
observed spectrum and (ii) the corresponding transmission spectrum through 
the atmospheric gases alone. The latter is calculated with a line-by-line RT 
model fed by the vertical compositional and thermal profiles predicted by the 
European Mars Climate Database for the dates, locations, and altitudes of the 
observations. Thus we define a significant new parameter /3 that expresses the 
strength of the gas-aerosols coupling while directly depending on r^^^. Com- 
bining /3 and the radiance value i?^^^ within the 2 /im band, we evaluate t^^^ 
and by radiative transfer inversion and provided that the radiative prop- 

erties of the aerosols are known from previous studies and that an independent 
estimation of the scale height of the aerosols is available. One should note 
that practically /3 can only be estimated for surfaces spectrally dominated by 
minerals and/or H2O ice. 

The validation of the method was performed both with synthetic and real 
data. It shows that our method is reliable if two conditions are fulfilled: (i) 
the observation conditions provide an airmass that is large enough > 3.5 (ii) 
the aerosol are vertically well mixed in the atmosphere. A low lying aerosol 
layer is not detectable by our method because it does not have much influence 
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on the mean pathlength of the photons into the gas. By contrast, a detached 
layer high in the atmosphere superimposed on the standard profile has a strong 
influence on the coupling. Our method works even if the underlying surface is 
completely made of minerals, corresponding to a low contrast between surface 
and atmospheric dust, while being observed at a fixed geometry. This the first 
principal asset of our method. Minimizing the effect of the surface Lambertian 
assumption on the AOD retrieval is the second principal asset of our method. 

Experiments conducted on OMEGA nadir looking observations -as well as 
CRISM EPF- produce very satisfactory results. With OMEGA, we note a good 



coherency between our approach and the one of Vincendon et al. (20081. In 
addition the coherence between the scale height deduced from the analysis of 
the observation and the value that we use a priori in the inversion process is 
verified in most cases. The domain of airmass for which our method is reliable 
implies that it is intended for analyzing high latitude OMEGA observations 
with sufficiently high solar zenith angles. This constraint is somewhat lighten 
with CRISM EPF observations that imply a large range of emergence angles 
and thus values of airmass. Indeed we note very satisfactory r^f^ retrievals for 
solar incidence angle down to 33°. Finally we should note that our method 
was applied for the first time extensively on a series of OMEGA observations in 
order to map the atmospheric dust opacity at high southern latitudes from early 
to late spring of Martian Year 27. This study is presented in the companion 
paper. 
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